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We investigate energetic stability of vortices near the surface of a Bose-Einstein condensate. From 
an energy functional of a rotating Bose-Einstein condensate, written in terms of variables local to the 
surface, and a suitable trial wavefunction we calculate the energy of a moving vortex. The energetic 
stability of the vortex is investigated in terms of the rotation frequency of the confining potential. 
The critical frequency at which the vortices enter the condensate is calculated and compared with 
the experiment where a reasonable agreement is obtained. 
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I. INTRODUCTION 

Since the experimental achievement of Bose-Einstein condensation in confined dilute gases and the observation of 
vortex lattice s 111 1^ L a, La la I fl , l7 | . much interest has been devoted to the formation, stabilization, and dynamics of 
vortices [il H I sTTaTl o T fill TiaTl la. liH H H ITil H7L liH ] . Vortices have been created in Bose-Einstein condensates using 
two different techniques. The first is a phase imprinting technique |l4j| and the second is a rotating harmonic trap 
technique [lj, ta, III li| • In the latter technique, a rotating laser beam superposed on the magnetic trap creates an 
anisotropic confining potential that stirs the condensate. Energetically stable vortices are generated if the rotation 
frequency f2 of the harmonic trap is greater than a certain critical frequency f2 c fl6l IT^. I20L . 

As a mechanism for vortex nucleation, it was su gge sted that vortices are created from the excitation of low-energy 
surface modes of the Bose-Einstein condensate [T^ . l2lll23.l23ll25l |. For rotational frequencies Q < f2 c , moving vortices 
are created in the ultra-dilute reg ion of the condensate. Their effect is observed at the surface of the condensate 
only as a surface wave |l7l Il8l 126. 28]. Upon increasing the rotation frequency, the vortices approach the surface of 
the condensate. At the critical frequency Q c , the vortices enter the condensate, i.e., nucleate. This was supported 
by Penckwitt et. al j2(J with a simulation to the growth of the condensate from a rotating thermal cloud. It was 
also supported by Anglin [rH ITsf who showed, using a boundary layer approach, that the energy barrier to vortex 
penetration disappears at the Landau critical velocity for surface modes. 

This picture for vortex nucleation will be the main focus of the present paper. Using a (rather simple) variational 
approach to the problem, we investigate the energetic stability of a vortex moving near the surface of the condensate. 
We verify the above picture by calculating the energy in terms of the distance of the vortex core from the surface Xq, 
the rotation frequency f2, and the center-of-mass speed of the core of the vortex vq. We also calculate the critical 
frequency f2 c at which the vortex enters the condensate. 

In terms of variables local to the surface of large condensates, the surface can be approximated by a plane and 
the harmonic trapping potential can be approximated by a linear potential [23, UK . We start by writing the Gross- 
Pitaevskii energy functional of a rotating condensate in this planar geometry. Using this energy functional with an 
appropriate trial wavefunction we calculate the energy of a vortex moving in a direction parallel to the surface of 
the condensate. Then we minimize the energy with respect to the position of the core of the vortex. The resulting 
equilibrium position of the vortex core is a function of the rotation frequency, the speed of the vortex, and the distance 
between the vortices in the case when more than a single vortex exists. This allows for a detailed investigation of the 
stability of the vortex in terms of these parameters. 

The rest of the paper is organized as follows. In Sec.|n]we introduce the energy functional of a rotating condensate 
in the planar geometry. In Sec. lIIII we use this energy functional to calculate the energy of a moving vortex as described 
above. In Sec. UVI we investigate the stability of the vortex in terms of the parameters involved. We end in Sec. IVIbv 
a summary and discussion of our results. 



II. ENERGY FUNCTIONAL OF THE SURFACE OF A ROTATING CONDENSATE 



In this section we write the Gross-Pitaevskii equation of a rotating condensate in terms of parameters local to the 
surface of the condensate. Then we construct an energy functional that corresponds to this equation. 
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The time-dependent Gross-Pitaevskii equation, that describes the behavior of the order parameter ip(r, t), is given 



2m 



V 2 + V(r,t)+g\^(r,t)\ 2 -p 



ip(r,t) = ih—tp(r,t) . 



(1) 



Here V(r,t) is a time-dependent harmonic trapping potential with time-dependence corresponding to rotating the 
axes of the trapping potential. The effective two-particle interaction g is proportional to the s-wave scattering length 
a according to g — AitaTi 2 /m, where m is the mass of an atom, and /i is the equilibrium chemical potential. 

We consider the trap axes in the x- and y-directions to be rotated counter clockwise around the z-axis with angular 
frequency f2. The time-dependence of the trapping potential can be removed by writing the Gross-Pitaevskii equation 
in the rotating frame 



-^-V 2 + V(T)+g\iP(r,t)\ 2 -n.L-» 
Zm 



d 

ip(r,t) = ih—ifj(r,t) 



(2) 



where L is the angular momentum operator. The trapping potential, which for simplicity is taken to be isotropic, is 
given by 



V(r) = -mw 2 r 2 , 

where luq is the characteristic frequency of the trap. Noting that f2L = —ihtt ■ (r x V) 
Gross-Pitaevskii equation takes the form 
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V 2 + V(r) + flr|V(r, t)\ 2 + i%(Sl x r) • V - p 



ip{v.t) = in—ip(Y,t) 



(3) 

-ih(n x r) ■ V, the 
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Next, we express the Gross-Pitaevskii equation in terms of the surface local coordinate r' and surface local velocity 
v' via a position transformation 



r = Rr + r' , 



(5) 



where R — -y/2/i/mcjQ is the Thomas-Fermi radius of the condensate. The Gross-Pitaevskii equation takes the form 



2m 



V 2 , + V(Rr + r') + g\ifj{r' , t)\ 2 + ih(fl x r') ■ V r - + ifi(n x Rr) ■ V r 



^(r , ,t)=ta^V(r',t). (6) 
at 



The term ih(Sl x Rr) ■ V r ' can be eliminated by a Galilean transformation to a frame of reference that is moving with 
a speed SI x RR, |29|. Thus, the Gross-Pitaevskii equation becomes 
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V 2 , + V(JZ? + r') + g\iP(r', t)\ 2 + ifi{VL x r') ■ V r 
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i>{r',t)=m—i){r',t) 
at 



(7) 



For distances close to the surface of the condensate, such that r' <C i?, the trapping potential can be approximated 
by a linear function as [2lJ 



V(Rr + r') ~ V{R) + V r ,V{Rr + r')\ r , =0 -r' 



V{R) + mu 2 R 
V(R) + Fx 



2ns J 



(8) 



where we have defined the force constant F — ttiuJqR and the surface local coordinate x = f • r', which is normal to 
the surface of the condensate. The coordinates of the y and z-components of r' are therefore parallel to the surface. 
In this plane geometry the Gross-Pitaevskii equation takes the form 



— V 2 , + Fx + g\^(r', t)\ 2 + ih((l x r') • V r 
2 m 



i)(r',t)=ih^(r',t) . 



(9) 



The first and the last terms on the left hand side of this equation can be combined as a complete square, which results 
in the following form for the Gross-Pitaevskii equation 



2m 



fl x r'J +Fx + g\ip(r',t)\ 
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m(n x r ) 



i)(r',t)=ih—i}(r',t) 



(10) 
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This equation can be derived by differentiating the following lagrangian with respect to ijj*(r' , t) 
where E[ip, ip*] is given by 

h 2 



= / dr' 
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V r ,-— fixr')^(r',i) 



+ Fx\^(r',t)\ 2 + i 5 |^(r',t)| 4 - \m{n x r') 2 |^(r', i)| 2 



(11) 
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III. ENERGY OF A MOVING VORTEX NEAR THE SURFACE OF THE CONDENSATE 

We consider a vortex with axis parallel to the surface of the condensate, say parallel to the z-direction. The core 
of the vortex is located at a distance x — xq from the Thomas- Fermi surface of the condensate and is moving in the 
positive y-direction with a speed Vq. This is depicted in the schematic figure Fig. ^ For simplicity we consider a 
vortex with a small core size such that the equilibrium condensate density can be considered constant over distances 
of the order of the core size. 



vortex core 




condensate surface 



FIG. 1: Figure showing the vortex core located at a distance Xo from the Thomas-Fermi surface of the condensate and is 
moving with a speed vo- 



The wavefunction of the condensate that contains such a vortex can thus be written as the product of the wave 
functions of the condensate that contains no vortex and the wave function of the vortex with an appropriate phase 
factor. The phase factor should take into account both the superfluid and the center of mass velocities of the vortex. 
Specifically, we write the wave function of the vortex as 



ip{r', t) = y / n (x)x(x -x ,y- v t) exp - x , y - v t)] 



(13) 



Here, Uq(x) = \ipo(x)\ 2 is the density profile of the vortex-free condensate at equilibrium. The equilibrium wavefunction 
i/jq(x) obeys Eq. ljTU|) after setting fl and dip/dt to zero, namely 



2m dx 2 



Fx + g\Mx)\ 2 



Mx) = . 



(14) 



The normalized density profile of the vortex is given by x 2 ( x ~ x o, V — Vot), where we assume for x(% — xq, y — Vq£) 
the trial function 30] 



x(x -x ,y- v t) 



(15) 



whith p = y/ {x — Xo) 2 + (y — Vot) 2 , and £ = 1/ ^/8Tran (x ) is the coherence length of the condensate calculated at 

Xq. 
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The phase of the vortex $(x — Xq, y — vot) must satisfy 



&tt * h 2 
— V r ^ = — + 

to mp 

= cos cp x 4- [Vo sin <p)y , 

mp mp 



(16) 



where p and are the polar coordinates in the x-y plane with the origin being at the core of the vortex as shown 
in Fig. ^ This equation expresses the fact that the velocity has two components: the superfluid velocity %<j)/mp 
associated with the vortex, and the center-of-mass velocity v^y. 

In Appendix A, we show that using the trial function Eq. <|13|) in the energy functional Eq. i|12|) . the energy of the 
vortex per unit length along the axis of the vortex takes the form 



E -- n (x ) / dp 
1 



r h 2 



x 



i 



n (x ) 



-mv. 



2m V P J 2 

HQ — m£l xqVq + F xq 



gn {xo)x 4 
dpx 2 



(17) 



The quantity of interest is the extra energy associated with the presence of the vortex |3l| . Therefore, we have to 
subtract the vortex energy given by the last equation from the energy of a vortex-free system that contains the same 
number of atoms as the system with the vortex. The energy density of the vortex-free system is only that of the 
mean-held interaction, the trapping potential, and the energy associated with the rotation, namely 



1 



gnl(x ) + F x n (x ) + Ti£ln (x ) , 



(18) 



where the subscript "vf" in E v { denotes vortex free. To fulfil the requirement that the number of atoms v in both 
systems must be the same, we calculate the energy of both systems within a cylindrical region of axis normal to the 
x-y plane and circular cross-section of radius b. The radius b must be, on the one hand, larger than several coherence 
lengths £ in order to incorporate all density variations associated with the presence of the vortex. On the other 
hand, b must be much smaller than the length over which the density changes significantly in order not to violate the 
local-density approximation used to derive calculate the energies Eqs. i|17|) and (|18|) . The number of atoms per unit 
length of the axis of the cylinder in the system that contains the vortex can be expressed as 



v = n (x ) / 2irpdpx 



Trb 2 n (xa) - n {x ) / 2-Kpdp(l~x 2 ) 



(19) 



Using this number to calculate the average density as no(xo) = v 2 /irb 2 , the energy per unit length of the vortex-free 
system is given by 



E v{ ~ ^tt b 2 gnl(xo) - gnl(x ) J 2npdp(l - x 2 ) + (Fx n (x ) + hnn (xo)) ^2i:pdpx 2 - (20) 

In this equation we ignore terms that are second order in 1 — x 2 hi the mean-field energy |3lj ] . 

Subtracting Eq. (J20J from Eq. JHJ, we obtain the energy associated with the presence of the vortex per unit length 
of its axis: 



E v = n (x ) J dp 



2m 



^[(Vx) 2 + h)+^gn (xo)(x 2 -l) 



1 



n (a;o) 



1 



-mv — mVl xqVq 



j dp X 2 + hflj dp( X 2 - 1) + Fx n (x ) J dp( X 2 - 1) 



(21) 



We recognize the first line of this equation as the internal energy of the vortex |31 
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IV. STABILITY ANALYSIS 



In this section we consider an infinite string of equidistant moving vortices. The string is parallel to the surface of 
the condensate, say along the y-axis, and the vortices move in that direction with a constant speed Vq. The distance 
between two adjacent vortices A is large enough such that the cores do not overlap. In this case, the energy of the 
string of vortices per unit length of the string is just the energy of one vortex given by Eq. I|21|) times the number of 
vortices per unit length. Therefore, it is sufficient to perform the analysis by calculating the energy of a single vortex. 
This energy will be a function of A that corresponds to a certain number of vortices per unit length. 

In the following, we investigate the dependence of the equilibrium position of the vortex on the rotation frequency 
fi, the distance between two adjacent vortices A, and the vortex speed vq. 

For convenience, we rewrite Eq. Ij21(l in terms of scaled quantities as follows: 



E v = e v + vin (x )(v$ - 2ttx v ) + v 2 n (x ){2tt + x ) , 



(22) 



where v\ — 2ir J Q b pdpx 2 , v 2 = 2ir J° pdp(x 2 ~ 1), and 



b _ 



e v = 27rn (So) / pdp 
Jo 



(Vx) 2 + |f 
P 



tUo(£o)(x 2 



(23) 



is the internal energy of the vortex. Here, the gradient operator is with respect to p. Length is scaled to the 
characteristic length at the surface of the condensate 5 defined by h 2 /2mS 2 = FS. The vortex velocity vo is scaled to 
U/mS, the angular frequency fl is scaled to h/mS 2 , the density n is scaled to 87ra<5 2 , and the energies E v and e v are 
scaled to (h 2 /2mS 2 ) / \8ira/5) — ((2i?/ao) 1 / 3 /167r)?icJo- In these dimensionless variables, Eq. 1(1^1 . takes the form 



dx 2 



n (i ) 



\/no(io) = . 



Furthermore, the vortex profile x, defined in Eq. (|15|) . takes the form 



(24) 



y/2/h (x ) + P 2 



(25) 



Since the distance between two adjacent vortices is A, the upper limit of the integrations in the expressions for v 2 , 
and Eq. §L\\\ is taken as b = A/2. 

We start by minimizing the vortex energy E v with respect to the vortex linear speed vq. This gives simply vq = xqQ, 
which is similar to the relation between the center of mass speed and the angular frequency of a rolling rigid wheel of 
radius xq. Substituting this value for vq in Eq. (|22(l . the energy of the vortex becomes 



EL 



viho(xa)x SI 2 + v 2 h {x ){2£l + x a ) 



(26) 



The equilibrium position of the vortex i oq is obtained by minimizing the energy E v with respect to the vortex 
distance from the surface xq, for a given set of the parameters Q and A. The main focus will be on the behavior of the 
equilibrium position a; eq as a function of the parameters f2 and A. In particular, we are interested in situations when 
the vortex enters the condensate. The Thomas-Fermi surface, defined by xq = 0, is to be taken as the boundary for 
considering the vortices inside or outside the condensate. Vortices located in the region xo > are to be considered 
out of the condensate, while those in the region xq < are considered inside the condensate. 

To calculate the vortex energy E v , we first need to solve Eq. (|24l) for n (io). As a first approximation we use 
the Thomas-Fermi approximation for the condensate density no(%o)- This is obtained by neglecting the first term in 
Eq. namely 



ho(x ) = 



-x , 
0, 



Xq < 
Xq > 



(27) 



The Thomas-Fermi approximation is accurate only in the deep region of the condensate xo <^C 0. Thus, the results 
we obtain here using this approximation will be accurate only in that region. From Eq. (|25|1 we see that the core size 
of the vortex depends on the condensate density. This will complicate the dependence of the energy on the vortex 
location. For simplicity we take the size of the core to be constant and is equal to its value at the Thomas-Fermi 
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surface. Substituting for ho(xo) from Eq. I|27[l in Eq. I|26[l . the vortex energy E v becomes a cubic function in Xq, 
which can be written as 

E v = -(«! + 2Qu 2 )x + (a 2 - v 2 )x\ + v x Vl 2 x% (28) 

where 



P 2 



and 

V2=2kJ o pd~p\{x 2 -I) 2 ■ (30) 

To be able to extract the main features of this energy expression, we point out some properties of the numbers 
v\ and v 2 and the parameters a.\ and a 2 . The number v\ is always positive since it corresponds to the number of 
atoms within a circle of radius b centered at the axis of a vortex that is embedded in a uniform background of density 
ho = 1. The number v 2 is always negative since it is the difference between the number of atoms in the vortex V\ 
and the number of atoms in a vortex-free background. The parameter a\ corresponds to the kinetic energy due to 
the gradient of the wavefunction of the vortex. The parameter a 2 corresponds to the mean-field interaction energy. 
It should be also noted that Eq. (|28|l is defined only for negative values of £q- 

With this in mind, we conclude from Eq. H28|) that E v starts linearly for small xq with a slope — (ai + 2Q,v 2 ) that, 
depending on the value of tt, can be either positive or negative. For large xo the energy E v behaves as Vi(1 2 Xq which 
is large and negative. If the initial slope — (a± + 2Qv 2 ) is positive, E v will have only a local maximum at some value 
of xq. If, on the other hand, the slope is negative, E v may have a local minimum in addition to a local maximum. If 
the negative value of the slope is large enough, the value of E v at the local maximum can vanish or even disappear. 
To be more explicit, we rewrite Eq. (|28|l as 

E v = x [(x - xi)(x - x 2 )} , (31) 

where 



•'•i ~~7~nr~ 1 ^To7\ M<*(.«i -:■ lV,<>). CV2) 



and 



*2 = ° 2 ~? ^rji^ v 2 ) 2 + 4^( Ql + 2v2 &) , (33) 

are two of the three roots of Eq. 131|) . The value of E v at the local maximum vanishes when X\ = x 2 . This can be 
satisfied if 

{a 2 -v 2 f +±v 1 n 2 {a l + 2v 2 n) = 0. i (34) 

which is the condition that gives the critical frequency f2 c at which the height of the energy barrier for vortex nuclcation 
vanishes. 

The expression for E v is plotted in Fig.|21as a function of Xq for different values of Q and one value of A. In addition 
to the above-mentioned properties of the energy curve, the figure shows that the height of the barrier decreases with 
increasing 17. At a certain value of fl the height of the barrier vanishes and the vortex can enter the condensate. For 
the value of A used in this figure, we find that, at f2 = 0.5, the height of the barrier vanishes. This critical value of Cl 
is indeed equal to the real root of Eq. H34JI ■ 

Since the Thomas-Fermi approximation, used here, is not accurate at the surface, we will not attempt to calculate 
the critical frequency at which the vortex enters the condensate. This is to be investigated next where we use the 
exact density profile of the condensate. 

For a more accurate treatment of the problem, we repeat the above calculation using the numerical solution of 
Eq. (|24|l instead of the Thomas- Fermi one. We start by solving numerically Eq. (|24fl for no(xo), and then using this 
result in Eq. 125|) to calculate Xi where the coherence length £ depends on ho(xo) and is not taken to be constant as in 
the above calculation. These values of fio(^o) and X are then used to calculate E v from Eq. (|26|l . Similar to the above 
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FIG. 2: Vortex energy E v as a function of the distance between the core of the vortex and the Thomas-Fermi surface xq. The 
energy is calculated using the Thomas-Fermi approximation. Five curves are plotted for values of Q that equal, starting from 
the uppermost curve, 0.1, 0.45, 0.5, 0.55, 0.6. The value of A used is 21. The unit of energy is ((2i?/ao) 1//3 /167r)?kJo and the 
unit of distance is 8. For a typical experiment such as the one discussed at the end of Sec. IV, the unit of energy equals 
0.033fttJo, and the unit of distance equals 0.6ao, where ao = ^Ti/moJa. 





FIG. 3: Vortex energy E v as a function of the distance between the core of the vortex and the Thomas- Fermi surface xq. The 
energy is calculated numerically. Five curves are plotted for values of Cl that equal, starting from the uppermost curve, 0.0, 
0.17, 0.34, 0.51, 0.68. The value of A used is 21. Note that the height of the energy barrier vanishes for = 0.34 



result obtained using the Thomas-Fermi approximation, we find that E v has a local minimum and a local maximum, 
and tends to — oo as Xq — > — oo. In Fig. El we plot the energy E v as a function of Xq for different values of f2. This 
figure is to be compared with Fig. 4 of Ref. 01 an d Fig. 5 of Ref. [l^. The main feature of Figs. [3 and [31 which 
is that an energy barrier for vortex nucleation exists and that the height of the barrier decreases with increasing Q, 
agrees with the above-mentioned two figures. However, a detailed comparison is not possible for the following reasons. 
In Fig. 4 of Ref. the energy is calculated in the outer region of the condensate, which means that it is accurate 
only for large distances away from the Thomas-Fermi surface. In Fig. 5 of Ref. the energy is calculated using the 
Thomas-Fermi approximation and therefore is accurate only inside the condensate. The present calculation is valid 
in the region near the surface of the condensate from both sides. Thus, Fig. |31 presents a calculation of the energy in 
a region where both calculations of Refs. 0] and p| are not applicable. 

In Fig. 21 we plot the position of the local minimum of the energy and the height of the energy barrier as a function 
of f2. Several curves are plotted corresponding to different values of A. Fig. 0Ja) shows that the vortex approaches 
the surface when Cl is increased. To enter the condensate, the energy barrier has to vanish. Fig. Efb) shows that the 
energy barrier does indeed vanish for a certain value of (l. This value is to be considered as the critical value f2 c at 
which the vortex nucleates. 

Now, we can compare our predictions for the critical frequency at which the vortex nucleates with those observed 
experimentally. To that end, we need to express Q c in terms of the harmonic trap frequency ujo for a spherical 
condensate. This is readily obtained by noticing that the characteristic length scale at the surface of the condensate 
6 is related to the Thomas-Fermi radius of the condensate through 6 — (fr 2 /2mF) 1 / 3 = (ao/2i?) 1 / 3 ao, where ao = 
y/h/majQ. Using this relation, can be expressed in terms of ui as Cl = (Q/ujo)(a /2R) 2 / 3 , which simply leads to 
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FIG. 4: a) The equilibrium position of the vortex core versus the rotation frequency of the magnetic trap. These curves are 
calculated for values of A that equal, starting from the lowest curve, 10, 12, 14, 16, 18, 20. The sudden drop in two of the 
curves appear because the local minimum of E v is no longer present. The local minimum and local maximum merge to form a 
stationary point, b) The height of the energy barrier calculated for values of A that equal, starting from the rightmost curve, 10, 
12, 14, 16, 18, 20. The units of energy and distance are as those mentioned in Fig. 2, and the unit of frequency is (2R/ao) 2 ^ 3 uio, 
which for the experiment discussed at the end of Sec. IV equals 2.6u>o- 



Q = {2R/a a ) 2/3 Quja. 

In the experiment of the ENS group a single vortex was created in about 1 x 10 5 atoms of a 87 Rb condensate 
at a rotating frequency — 2ir x 147 Hz. The magnetic trap used was axially symmetric with axial trap frequency 
u> z = 2ir x 11.7 Hz and radial trap frequency u>± = 2ir x 219 Hz. The scattering length is taken as a — 5.8nm. For 
these parameters R ~ 2.1ao and (2i?/ao) 2//3 ~ 2.6, where we have used a>o = (h/m(uj'j_uj z ) 1 / 3 ) 1 / 2 . For a single vortex 
the value of b is half the circumference of the condensate, namely b « 2.l7rao « 10.5(5. Since A = 26, the lowest curve 
of Fig. 01a) and the left most curve of Fig. 0^b) correspond to the case of a single vortex. For the latter curve the 
height of the barrier vanishes at fl c ss 0.34. (Notice that this corresponds to the number obtained above using the 
Thomas- Fermi approximation, namely Q c ~ 0.5). Thus, according to our calculation, the critical frequency for vortex 
nucleation is f2 c /^o = 0.88. This is to be compared with the experimental number H c /uj± = 147/219 = 0.68. Notice 
that we consider c^o to correspond to u± since in the experiment the axis of rotation is the z-axis. Taking into account 
the approximations that we have used, we consider this agreement to be reasonable. One possible reason for the 
disagreement is that we have taken the core of the vortex to be small compared to the distance over which the density 
of the condensate changes significantly. This approximation is accurate only in the deep region of the condensate and 
not close to the Thomas-Fermi surface or away from it. Furthermore, our criterion for vortex nucleation may be too 
restrictive since, even with a nonzero barrier height, as shown in Fig. E^b), the vortex can tunnel through the barrier 
and thus nucleate at a frequency which is less than that obtained above. 



V. CONCLUSION 



We have written down the energy functional of a rotating Bose-Einstein condensate in terms of variables local to 
its surface. Using a variational ansatz, we have calculated the energy of a single vortex moving at the surface of 
the condensate. This energy was used to investigate the energetic stability of the vortex in terms of the rotation 
frequency of the confining potential. Upon increasing the rotation frequency, the vortex approaches the surface of the 
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condensate. At a critical rotation frequency, the vortex crosses the surface of the condensate and enters it. 

The local-density approximation was used to simplify the calculations and lead to analytical results in the Thomas- 
Fermi limit. This approximation is accurate only in the interior of the condensate and not near the Thomas- Fermi 
surface or away from it. The critical frequency obtained using this approximation is thus not expected to be accurate. 
Therefore, we solve the problem numerically including kinetic energy effects at the surface. This leads to more accurate 
values for the critical frequency that compare reasonably with experiment. It should be also noted that our criterion 
for vortex nucleation, which is requiring that the energy barrier vanish, may be too restrictive. When the height of the 
energy barrier is nonzero but small enough, the vortex can tunnel through it and thus nucleate at a rotation frequency 
that is less than the one we calculated in the previous section. In this sense, our calculated critical frequency may be 
regarded as an upper bound to the vortex nucleation frequency. 



APPENDIX A: DETAILS OF THE ENERGY CALCULATION 



In this Appendix we present details for the derivation of Eq. I|17|) from the energy functional Eq. (|12|l using the 
trial wave function Eq. i|13fl for a moving vortex. 

Since the vortex axis is parallel to the z-axis and the vortex is moving in the y-direction, the nontrivial contributions 
to the energy arise only from integrating over x and y. Therefore, we calculate the energy of the vortex per unit 
length of the z-direction, which means that we have to perform only integrals in the x-y plane. It is, thus convenient 
to shift the origin of r' from the surface of the condensate to the core of the vortex and to express the integral of 
energy functional in the polar coordinates p and (f>. Explicitly, we make the transformation r' = xqH + p. 

The kinetic energy contribution reads 



Ek = / dp— 

2m 



V - — SI x (x x + p) ) ip(p, t 



(Al) 



where, here and through out this appendix, the gradient term is with respect to p. Inserting the trial function Eq. i|13|) 
in this equation we obtain 



E k - ^ n o( x o) I dp 



V - — SI x (x x + p) ) X (p) exp (i$(p)) 



(A2) 



Here, we have employed the approximation that the condensate equilibrium density no(x) is constant over distances of 
the order of the coherence length. For this approximation to be valid, the upper limit of the integration over p must be 
much less than the distance over which the condensate density changes significantly. In this case we can approximate 
no(x) by its value at xq and take it out of the integration. Furthermore, the last equation can be simplified to 



h 2 f 
Ek = 2^ n °^ J df> 



V X + ix V $ — Si x (x x + p)x 



(V%) 2 + (xV$)^ + ^ (fl x (x x + p) f x 2 - 2^ (fi x (x x + p)) ■ V$ 



(A3) 



Knowing that O = flz and using the expression for 4> from Eq. lltjfl . it is straightforward to show that the kinetic 
energy takes the form 



E k 



2m 



n {x ) / dp 



(Vx) s 



x_ 

P 2 



2„,2 



+ 



m v. 



o ,,2 



2m 



X 



x v )x 2 + jTn 2 (x 2 +p 2 ) X 2 



The contribution to the energy from the last term of Eq. I|12fl is 

1 

En = 



■-mno(xo) / dp (SI x (x x + p)) x 
= -\mSl 2 71q(xq) I dp(xQ + p 2 )x 



(A4) 



(A5) 



When added to the kinetic energy, this term cancels with the last term of Eq. (|A4|) . 
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The contribution of the external trapping potential is given by 

E oxt = Fn (x ) J dp(x + /9sin0)x 2 

= Fn (x )x / dp X 2 ■ (A6) 



The contribution of the interatomic interaction reads simply 

E in t = \gn 2 (x ) Jdp X A . (A7) 
Finally, we add all these contributions to obtain Eq. (|T7Jl . 
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